#ifndef _DIRECTPOSSIONSOLVER_H_
#define _DIRECTPOSSIONSOLVER_H_
#include "EmbeddedGrid.H"
#include "FuncFiller.H"
#include <string.h>
#include <fstream>
#include <lapacke.h>

template <int Dim>
class DirectPossionSolver;

template <>
class DirectPossionSolver<2>
{
public:
  // Some convenient definitions
  enum
  {
    Dim = 2
  };
  using iVec = Vec<int, Dim>;
  using rVec = Vec<Real, Dim>;
  using EmGrid = ExcludedCircle<Dim>;
  template <class T>
  using vector = std::vector<T>;

public:
  /**
   * @brief Construct a new Direct Possion Solver object
   * @param  aGrid            My Param doc
   * @param  aBcTypes         My Param doc
   */
  DirectPossionSolver(const EmGrid &aGrid,
                      const char *aBcTypes,
                      const int ordIntpDiriIrBdr,
                      const int ordIntpNeumIrBdr);
  template <class TFunc>
  Tensor<Real, 1> computeSolution(const Tensor<Real, Dim - 1> *aRectBdrData,
                                  const Tensor<Real, Dim> &aDomainData,
                                  const TFunc &aIrBdrDataExpr) const;

  template <class VecExpr, class ScalerExpr>
  Tensor<Real, 1> computeSolution(const ScalerExpr &negDDexpr,
                                  const VecExpr &Dexpr,
                                  const ScalerExpr &expr) const;

protected:
  EmGrid emGrid;
  char bcTypes[5];
  int numUnknownsInDomain;
  int numUnknownsExceptIrGhost;
  int numUnknownsAll;
  Tensor<int, 2> unknownMap;
  vector<iVec> unknownMapInv;
  Tensor<Real, 2> poissonMat;

  // 控制不规则边界的精度
  //目前支持 1, 2
  int ordIntpDiriIrBdr;
  //目前支持 1, 2
  int ordIntpNeumIrBdr;

  void doFindUnknownMap(); // Step 1
  void doFillMat();        // Step 2
  void doFillMat_DiriIrBdr(const iVec &idx);
  void doFillMat_NeumIrBdr(const iVec &idx);

  template <class TFunc>
  const Tensor<Real, 1> doFillRhs(const Tensor<Real, Dim - 1> *aRectBdrData, const Tensor<Real, Dim> &aDomainData, const TFunc &aIrBdrDataExpr) const;
  const Tensor<Real, 1> doDirectSolve(const Tensor<Real, Dim - 1> &rhs) const;

public:
  /**
   * @param  q       q-norm of grid function. max-norm when q=0.
   */
  template <class TFunc>
  const Real computeErrNorm(const Tensor<Real, 1> &sol, const TFunc &exactSolution, const int q = 1) const;
  /**
   * @brief 将计算域网格上的数值解,真解与误差以方阵形式打印到指定文件中。
   * 空缺部分用0占位。
   * @param  _file  给定文件名。          My Param doc
   */
  template <class TFunc>
  const void plotIn(const Tensor<Real, 1> &sol, const std::string &_file, const TFunc &exactSolution) const;
};

DirectPossionSolver<2>::DirectPossionSolver(const EmGrid &aGrid,
                                            const char *aBcTypes,
                                            const int _ordIntpDiriIrBdr,
                                            const int _ordIntpNeumIrBdr)
    : emGrid(aGrid), ordIntpDiriIrBdr(_ordIntpDiriIrBdr), ordIntpNeumIrBdr(_ordIntpNeumIrBdr)
{
  strncpy(bcTypes, aBcTypes, 5);
  const auto dx = emGrid.spacing();
  // Step 1. Find 1-to-1 correspondence between mesh points and unknowns of matrix.
  doFindUnknownMap();
  // Step 2. Fill in the coefficients of matrix.
  doFillMat();
}

template <class TFunc>
Tensor<Real, 1> DirectPossionSolver<2>::computeSolution(const Tensor<Real, Dim - 1> *aRectBdrData,
                                                        const Tensor<Real, Dim> &aDomainData,
                                                        const TFunc &aIrBdrDataExpr) const
{
  // Step 1. Fill in the coefficients of Rhs.
  const Tensor<Real, 1> rhs = doFillRhs(aRectBdrData, aDomainData, aIrBdrDataExpr);
  // Step 2. Use lapacke to solve the equation.
  return doDirectSolve(rhs);
}

template <class VecExpr, class ScalerExpr>
Tensor<Real, 1> DirectPossionSolver<2>::computeSolution(const ScalerExpr &negDDexpr,
                                                        const VecExpr &Dexpr,
                                                        const ScalerExpr &expr) const
{
  FuncFiller<Dim> ff(emGrid);
  Box<Dim> gbx = emGrid;
  Tensor<Real, Dim> phi(gbx);
  Tensor<Real, Dim - 1> bcData[2 * Dim];
  loop_box_2(gbx, i, j)
      ff.fillDomainVal(phi, negDDexpr);
  for (int d = 0; d < Dim; ++d)
  {
    for (int side : {-1, 1})
    {
      int m = 2 * d + (side + 1) / 2;
      bcData[m].resize(reduce(gbx, m / 2));
      switch (bcTypes[m])
      {
      case 'D':
      {
        ff.fillRectBdrVal(bcData[m], d, side, expr);
        break;
      }
      case 'N':
      {
        ff.fillRectBdrNormalGrad(bcData[m], d, side, Dexpr);
        break;
      }
      default:
        break;
      }
    }
  }
  switch (bcTypes[4])
  {
  case 'D':
  case '\0':
    return computeSolution(bcData, phi, expr);
    break;
  case 'N':
    return computeSolution(bcData, phi, [&](const auto &x)
                           { return dot(Dexpr(x), (*this).emGrid.embeddingCircle.getNormalGrad(x)); });
    break;
  default:
    assert(!"To be implemented!");
    return (Tensor<Real, 1>)0;
    break;
  }
}

#endif